module test_mod
    implicit none
    
contains
    function func_test(e,nu,t) result(dm)
        real, intent(in) :: e,nu
        real, intent(inout) :: t
        real :: a
        real, dimension(3,3) :: dm
        t = t + 100
        a=e/(1-nu*nu)
        dm(1,1)=a
        dm(2,2)=a
        dm(1,2)=nu*a
        dm(2,1)=nu*a
        dm(3,3)=0.5*(1-nu)*a
    end function func_test
end module test_mod

program test_
    ! call test14()
    ! integer::n
    ! real,allocatable::a(:)
    ! n=10
    ! call test10(n,a)
    ! print *, "a =",a
    ! do n=1,10
    !     print *, "a(",n,")=",a(n)
    ! end do
    call test17()
contains

subroutine test17()
    use iso_fortran_env,only:real64
    implicit none
    integer :: i,ios
    real(real64)::arr(3),earr(3)
    character(len=256) :: str
    earr=[1.2509373668496_real64,-4.0574914942688_real64,3.1951856799543_real64]
    str='    1.2509373668496         -4.0574914942688          3.1951856799543   '
    str=trim(adjustl(str))
    print * , 'str=',str,'++'
    print *, 'len(str)=',len(str)
    print *, 'after reading...'
    
    do i=1,3
        read(str, *, iostat=ios) arr  ! 读取第i个数值到数组
        if (ios /= 0) exit
        print '(a,f40.35)','error=',arr(i)-earr(i)
    end do
end subroutine test17

subroutine test16()
    implicit none
    integer, parameter :: N = 10
    integer, parameter :: k = 4
    integer, dimension(k) :: b = [2, 5,5, 7]
    integer, dimension(:), allocatable :: res
    logical, dimension(:), allocatable :: mask   ! <-- 关键改动
    integer :: i, cnt

    ! 动态分配
    allocate(mask(N))
    mask = .true.
    do i = 1, k
        if (b(i) >= 1 .and. b(i) <= N) mask(b(i)) = .false.
    end do

    cnt = count(mask)
    allocate(res(cnt))

    cnt = 0
    do i = 1, N
        if (mask(i)) then
        cnt = cnt + 1
        res(cnt) = i
        end if
    end do

    write(*, '(10I6)') res

    deallocate(mask, res)
end subroutine test16

subroutine test14()
    real::a(3,3),b(3),x(3)
    integer:: i,j
    do i=1,3
        do j=1,3
            a(i,j)=i+j
            end do
        b(i)=3*i
    end do  

    print *, "a ="
    write(*,'(3F5.1)') transpose(a)

    print *, "b ="
    write(*,'(3F5.1)') b
    b=matmul(a,b)
    print *, "x ="
    write(*,'(3F5.1)') b
    
    print *, "b(1:3) ="
    write(*,'(3F5.1)') b(1:3)

    print *, "100*matmul(a,b)="
    write(*,'(3F10.1)') 100*matmul(a,b)
end subroutine test14

! subroutine test13()
!     ! use writer_mod 
!     implicit none
!     integer :: i,j,elems(20,4)
!     real:: nodes(10,3),nu(10,6),reac(10,6)
!     do i=1,10
!         do j=1,3
!             nodes(i,j)=i*j
!         end do

!         do j=1,6
!             nu(i,j)=i-j
!             reac(i,j)=10*i+j
!         end do 

!         do j=1,4
!             elems(i,j)=i*j
!         end do
!     end do
!     call ls_writer(nodes,10,elems,20,nu,reac)

! end subroutine test13

subroutine test12()
    integer ::a(15,15),i,j
    do i=1,15
        do j=1,15
            a(j,i)=i
        end do
    end do

    print *, "a ="
    do i=1,15
        write(*,'(15I4)') a(i,:)
    end do

    print *, "a(1,1::3) ="
    print *, a(1,1::3)

end subroutine test12

subroutine test11()
    implicit none
    integer:: el_conn(4)=[11,24,34,4]
    integer:: el_inds(24)
    integer:: i,j
    do i=1,4
        j=el_conn(i)
        el_inds(6*i-5)=6*j-5
        el_inds(6*i-4)=6*j-4
        el_inds(6*i-3)=6*j-3
        el_inds(6*i-2)=6*j-2
        el_inds(6*i-1)=6*j-1
        el_inds(6*i)=6*j
    end do

    write(*,'(A,/,24I5)') "el_inds =", el_inds
end subroutine test11

subroutine test10(n,ou)
    integer,intent(in) :: n
    real,allocatable,intent(out) :: ou(:)
    integer::i

    allocate(ou(n))
    do i=1,n
        ou(i)=i
    end do
end subroutine test10

subroutine test9()
    real::k(3,3)
    integer::i,j
    call random_seed()
    do i=1,3
        do j=1,3
            call random_number(k(i,j))
        end do
    end do
    k=k*100
    print *, "k ="
    do i = 1,3
        write(*,'(3F10.4)') k(i,:)
    end do

    print *, "k(:,1) ="
    write(*,'(3F10.4)') k(:,2)

    print *, "size of k(:,1) =", shape(k(:,1))
    
end subroutine test9

! subroutine test8()
!     use constants_mod, only: kme_inds,kbs_inds
!     real::kme(12,12),ke(24,24)
!     integer::i,j
!     ke=0.0
!     call random_seed()
!     do j=1,12
!         do i=1,12
!             call random_number(kme(i,j))
!         end do
!     end do

!     do i=1,12
!         write(*,'(12F8.4)') kme(i,:)
!     end do

!     ke(kme_inds,kme_inds)=kme
!     do i=1,24
!         write(*,'(24F8.4)') ke(i,:)
!     end do
! end subroutine test8

subroutine test7()
    implicit none
    integer,parameter::n=5
    integer::i,j
    real:: a(n,n),b(n,n),c(n,n)
    character(len=15) :: fmt
    
    call random_seed()
    do j=1,n
        do i=1,n
            call random_number(a(i,j))
            call random_number(b(i,j))
            call random_number(c(i,j))
        end do
    end do
    write(fmt,'(A,I0,A)') '(',n,'F8.4)'
    print *, "fmt=",fmt
    print *, "origin array a:"
    do i=1,n
        write(*,fmt) a(i,:)
    end do
    print *, "origin array b:"
    do i=1,n
        write(*,fmt) b(i,:)
    end do
    print *, "origin array c:"
    do i=1,n
        write(*,fmt) c(i,:)
    end do

    print *, "after c=matmul(transpose(b),matmul(a,b)) :"
    c=matmul(transpose(b),matmul(a,b))
    do i=1,n
        write(*,fmt) c(i,:)
    end do

end subroutine test7

subroutine test6()
    integer :: i,j,k
    integer::arr(4,4),sub_arr(2,2),a42(4,2)
    k=0
    do i=1,4
        do j=1,4
            k=k+1
            arr(i,j)=k
        end do
    end do
    
    print *, "origin array:"
    do i=1,4
        write(*,'(4I4)') arr(i,:)
    end do
    
    sub_arr=arr(2:3,2:3)
    print *, "sub array, arr(2:3,2:3):"
    do i=1,2
        write(*,'(2I4)') sub_arr(i,:)
    end do

    sub_arr=arr([1,3],[1,3])
    print *, "sub array, arr([1,3],[1,3]) :"
    do i=1,2
        write(*,'(2I4)') sub_arr(i,:)
    end do

    a42=arr(:,[1,4])
    print *, "sub array, arr(:,[1,4]) :"
    do i=1,4
        write(*,'(2I4)') a42(i,:)
    end do

end subroutine test6

! subroutine test_inverse()
!     use matrix_utils_mod
!     implicit none
!     real :: A(3,3), A_inv(3,3),B(3,3)
!     integer :: i,ok
    
!     ! 测试矩阵
!     A = reshape([2.0, 0.00001, 0.0, &
!                     0.0, 3.0, 1555.0, &
!                     900.0, 0.0, 2.0], [3,3])
    
!     call invert_matrix(A, A_inv, ok)
    
!     if (ok==0) then
!         print *, "原始矩阵:"
!         do i = 1, 3
!             print '(3F10.4)', A(i,1),A(i,2),A(i,3)
!         end do
        
!         print *, "逆矩阵:"
!         do i = 1, 3
!             print '(3F10.4)', A_inv(i,:)
!         end do
        
!         ! 验证 A*A⁻¹ ≈ I
!         print *, "验证 A*A⁻¹:"
!         B=matmul(A, A_inv)
!         do i = 1, 3
!             print '(3F10.4)', B(i,:)
!         end do
!     else
!         print *, "矩阵求逆失败，矩阵可能奇异"
!     end if
! end subroutine test_inverse

subroutine test1()
    implicit none
    integer ::i, box(3,4)=reshape([1,2,3,4,5,6,7,8,9,10,11,12],shape(box))
    write(*,'(4I4)') transpose(box)

    do i=1,3
        write(*,'(4I4)') box(i,:)
    end do
end subroutine test1

subroutine test2()
use test_mod
    implicit none
    real :: e,nu,t
    real,dimension(3,3) :: dm
    integer :: arr(4,2),i,j
    dm(:,:)=0.0
    e=100000.0
    nu=0.3
    t=0.0
    arr=reshape([ 1,-1,-1, 1, 1, 1,-1,-1], [4, 2])
    ! arr=reshape([ 1,2,3, 4, 5, 6,7,8], [4, 2])
    dm=func_test(e,nu,t)
    do i=1,4
        write(*,*) (arr(i,j), j=1,2)
    end do
    print *, sizE(arr,1),size(arr,2)
    
end subroutine test2

subroutine test3()
    use, intrinsic :: ISO_C_BINDING
    implicit none
    integer, parameter :: n = 3, nrhs = 1
    real(8) ,target :: a(n,n), b(n),c(n,n)
    integer :: ipiv(n), info, i
    real(8), pointer :: ptr_a(:,:), ptr_c(:,:)
    type(c_ptr) :: a_loc, c_ol
    ! 系数矩阵
    a = reshape([2.0d0, 1.0d0, 1.0d0,3.0d0, 2.0d0, 1.0d0,2.0d0, 1.0d0, 2.0d0], [n,n],order=[1,2])
    
    ! 右端向量
    b = [4.0d0, 7.0d0, 5.0d0]
    write(*,*) " size b=",shape(b)
    
    print *, "before solve AX=B, a ="
    do i=1,n
        write(*,'(3F10.4)') a(i,:)
    end do

    c=a
    print *, "before solve AX=B, c ="
    do i=1,n
        write(*,'(3F10.4)') c(i,:)
    end do

    ! 解方程
    call dgesv(n, nrhs, c, n, ipiv, b, n, info)
    
    if (info == 0) then
        print *, "Solution:"
        do i = 1, n
            print '(F10.4)', b(i)
        end do
        
        print *, " after solving AX=B, a ="
        do i=1,n
            write(*,'(3F10.4)') a(i,:)
        end do

        print *, " after solving AX=B, c ="
        do i=1,n
            write(*,'(3F10.4)') c(i,:)
        end do
    else
        print *, "Error in DGESV, info =", info
    end if
    ptr_a=>a
    ptr_c=>c
    a_loc=c_loc(ptr_a)
    c_ol=c_loc(ptr_c)
    ! print *, "address of a =", a_loc
    ! print *, "address of c =", c_ol
end subroutine test3

subroutine test4()
    implicit none
    ! notice the use of "shape(box)" on the RHS
    integer :: box(3,4)=reshape([1,2,3,4,5,6,7,8,9,10,11,12],shape(box))
    ! integer,allocatable :: v(:,:)
    ! integer :: rc(2)
    ! basics0
    ! what is the current shape of the array?
    call printi('shape of box is ',box)
    ! change the shape
    call printi('reshaped ',reshape(box,[2,6]))
    call printi('reshaped ',reshape(box,[4,3]))

    !    ! fill in row column order using order
    !     v=reshape([1,2,3,4,10,20,30,40,100,200,300,400],[1,12])
    !     call printi('here is some data to shape',v)
    !     call printi('normally fills columns first ',reshape([v],[3,4]))
    !     call printi('fill rows first', reshape([v],[3,4],order=[2,1]))

    !     ! if we take the data and put in back in filling
    !     ! rows first instead of columns, and flipping the
    !     ! height and width of the box we not only fill in
    !     ! a vector using row-column order we actually
    !     ! transpose it.
    !     rc(2:1:-1)=shape(box)
    !     ! copy the data in changing column number fastest
    !     v=reshape(box,rc,order=[2,1])
    !     call printi('reshaped and reordered',v)
    !     ! of course we could have just done a transpose
    !     call printi('transposed',transpose(box))

    !    ! making the result bigger than source using pad
    !     v=reshape(box,rc*2,pad=[-1,-2,-3],order=[2,1])
    !     call printi('bigger and padded and reordered',v)
        
end subroutine test4

subroutine printi(title,arr)
    implicit none

    !@(#) print small 2d integer arrays in row-column format

    character(len=*),parameter :: all='(*(g0,1x))' ! a handy format
    character(len=*),intent(in)  :: title
    integer,intent(in)           :: arr(:,:)
    integer                      :: i
    character(len=:),allocatable :: biggest

    print all
    print all, trim(title),':(',shape(arr),')'  ! print title
    biggest='           '  ! make buffer to write integer into
    ! find how many characters to use for integers
    write(biggest,'(i0)')ceiling(log10(real(maxval(abs(arr)))))+2
    ! use this format to write a row
    biggest='(" > [",*(i'//trim(biggest)//':,","))'
    ! print one row of array at a time
    do i=1,size(arr,dim=1)
        write(*,fmt=biggest,advance='no')arr(i,:)
        write(*,'(" ]")')
    end do

end subroutine printi

end program test_

